We present new evolutionary synthesis models of M82 based mainly on observations consisting of near-infrared integral field spectroscopy and mid-infrared spectroscopy. The models incorporate stellar evolution, spectral synthesis, and photoionization modeling, and are optimized for 1-45 micron observations of starburst galaxies. The data allow us to model the starburst regions on scales as small as 25 pc. We investigate the initial mass function (IMF) of the stars and constrain quantitatively the spatial and temporal evolution of starburst activity in M82. We find a typical decay timescale for individual burst sites of a few million years. The data are consistent with the formation of very massive stars (> 50-100 Msun) and require a flattening of the starburst IMF below a few solar masses assuming a Salpeter slope at higher masses. Our results are well matched by a scenario in which the global starburst activity in M82 occurred in two successive episodes each lasting a few million years, peaking about 10 and 5 Myr ago. The first episode took place throughout the central regions of M82 and was particularly intense at the nucleus while the second episode occurred predominantly in a circumnuclear ring and along the stellar bar. We interpret this sequence as resulting from the gravitational interaction M82 and its neighbour M81, and subsequent bar-driven evolution. The short burst duration on all spatial scales indicates strong negative feedback effects of starburst activity, both locally and globally. Simple energetics considerations suggest the collective mechanical energy released by massive stars was able to rapidly inhibit star formation after the onset of each episode.